R Notebook

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✓ ggplot2 3.3.4     ✓ purrr   0.3.4
✓ tibble  3.1.2     ✓ dplyr   1.0.7
✓ tidyr   1.1.3     ✓ stringr 1.4.0
✓ readr   1.4.0     ✓ forcats 0.5.1
── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

library(MatrixGenerics)

Loading required package: matrixStats

Attaching package: ‘matrixStats’

The following object is masked from ‘package:dplyr’:

    count


Attaching package: ‘MatrixGenerics’

The following objects are masked from ‘package:matrixStats’:

    colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse, colCounts, colCummaxs, colCummins, colCumprods, colCumsums, colDiffs, colIQRDiffs, colIQRs,
    colLogSumExps, colMadDiffs, colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats, colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
    colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads, colWeightedMeans, colWeightedMedians, colWeightedSds, colWeightedVars, rowAlls, rowAnyNAs, rowAnys,
    rowAvgsPerColSet, rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods, rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps, rowMadDiffs, rowMads,
    rowMaxs, rowMeans2, rowMedians, rowMins, rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks, rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs,
    rowVars, rowWeightedMads, rowWeightedMeans, rowWeightedMedians, rowWeightedSds, rowWeightedVars

library(SingleCellExperiment)

Loading required package: SummarizedExperiment
Loading required package: GenomicRanges
Loading required package: stats4
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: ‘BiocGenerics’

The following objects are masked from ‘package:parallel’:

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ, clusterExport, clusterMap, parApply, parCapply, parLapply, parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from ‘package:dplyr’:

    combine, intersect, setdiff, union

The following objects are masked from ‘package:stats’:

    IQR, mad, sd, var, xtabs

The following objects are masked from ‘package:base’:

    anyDuplicated, append, as.data.frame, basename, cbind, colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep, grepl, intersect,
    is.unsorted, lapply, Map, mapply, match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
    table, tapply, union, unique, unsplit, which.max, which.min

Loading required package: S4Vectors

Attaching package: ‘S4Vectors’

The following objects are masked from ‘package:dplyr’:

    first, rename

The following object is masked from ‘package:tidyr’:

    expand

The following objects are masked from ‘package:base’:

    expand.grid, I, unname

Loading required package: IRanges

Attaching package: ‘IRanges’

The following objects are masked from ‘package:dplyr’:

    collapse, desc, slice

The following object is masked from ‘package:purrr’:

    reduce

Loading required package: GenomeInfoDb
Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with 'browseVignettes()'. To cite Bioconductor, see 'citation("Biobase")', and for packages 'citation("pkgname")'.


Attaching package: ‘Biobase’

The following object is masked from ‘package:MatrixGenerics’:

    rowMedians

The following objects are masked from ‘package:matrixStats’:

    anyMissing, rowMedians

lseq <- function(from, to, length.out){
  exp(seq(log(from), log(to), length.out = length.out))
}
theme_set(cowplot::theme_cowplot())

library(org.Hs.eg.db)

Loading required package: AnnotationDbi

Attaching package: ‘AnnotationDbi’

The following object is masked from ‘package:dplyr’:

    select

human_cell_cycle_genes <- select(org.Hs.eg.db, keytype = "GOALL", keys = "GO:0007049", columns = "ENSEMBL")[, "ENSEMBL"]

Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
'select()' returned 1:many mapping between keys and columns

library(org.Mm.eg.db)

mouse_cell_cycle_genes <- select(org.Mm.eg.db, keytype = "GOALL", keys = "GO:0007049", columns = "ENSEMBL")[, "ENSEMBL"]

'select()' returned 1:many mapping between keys and columns

mu_sup <- lseq(1e-4, 1e6, length.out = 1001)
poisson_pred <- cross_df(list(mu = mu_sup, factor = 10^seq(-8, 8))) %>%
  mutate(var = mu * factor)

gampoi_pred <- cross_df(list(mu = mu_sup, factor = 10^seq(-2, 2, by = 2))) %>%
  mutate(var = mu + mu^2 * factor) 

Prepare Data

Download data

# Work around for some bug in zellkonverter
# https://github.com/theislab/zellkonverter/issues/45
setAs("dgRMatrix", to = "dgCMatrix", function(from){
  as(as(from, "CsparseMatrix"), "dgCMatrix")
})

if(! file.exists("../data/klein_2015.h5ad")){
  download.file("https://data.caltech.edu/tindfiles/serve/f0d567c5-cea6-4a60-923e-e9fb4a4019e8/", "../data/klein_2015.h5ad")
}
if(! file.exists("../data/svensson_2017_1.h5ad")){
  download.file("https://data.caltech.edu/tindfiles/serve/3f89d3a5-6ceb-486e-95d4-9bd3f511a706/", "../data/svensson_2017_1.h5ad")
}
if(! file.exists("../data/svensson_2017_2.h5ad")){
  download.file("https://data.caltech.edu/tindfiles/serve/16dab9ea-4447-4e23-9aad-e68d052fd789/", "../data/svensson_2017_2.h5ad")
}
if(! file.exists("../data/nih3t3.h5ad")){
  download.file("https://data.caltech.edu/tindfiles/serve/a448e98e-89cd-47b3-a134-803bbde29781/", "../data/nih3t3.h5ad")
}
if(! file.exists("../data/hek293t.h5ad")){
  download.file("https://data.caltech.edu/tindfiles/serve/b2758046-9247-43ab-b8f0-68882b4f39a3/", "../data/hek293t.h5ad")
}
if(! file.exists("../data/nci-h1975.Rds")){
  download.file("https://github.com/LuyiTian/sc_mixology/raw/master/data/csv/sc_10x.metadata.csv.gz", "../data/nci-h1975-metadata.csv.gz") 
  meta <- read.csv("../data/nci-h1975-metadata.csv.gz")
  meta$cell_id <- rownames(meta)
  
  download.file("https://github.com/LuyiTian/sc_mixology/raw/master/data/csv/sc_10x.count.csv.gz", "../data/nci-h1975.csv.gz") 
  count_mat <- as.matrix(read.csv("../data/nci-h1975.csv.gz"))
  gene_info <- AnnotationHub::AnnotationHub()[["AH53537"]] %>%
    as.data.frame() %>%
    group_by(gene_id) %>%
    summarize(chromosome = dplyr::first(seqnames),
              gene_name = dplyr::first(gene_name),
              gene_biotype = dplyr::first(gene_biotype))
  row_df <- tibble(gene_id = rownames(count_mat)) %>%
    left_join((gene_info), by = "gene_id") %>%
    as.data.frame()
  saveRDS(SummarizedExperiment(S4Vectors::SimpleList(counts = count_mat), colData = meta, rowData = row_df), "../data/nci-h1975.Rds")
}
if(! file.exists( "../data/GSE126321.Rds")){
  download.file("https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE126321&format=file", "../data/GSE126321_RAW.tar")
  dir.create("../data/GSE126321")
  untar("../data/GSE126321_RAW.tar", exdir = "../data/GSE126321")
  mat <- Matrix::readMM("../data/GSE126321/GSM3596320_GM18502_matrix.mtx.gz")
  genes <- read_tsv("../data/GSE126321/GSM3596320_GM18502_genes.tsv.gz", col_names = c("gene_id", "gene_name"))
  barcodes <- read_tsv("../data/GSE126321/GSM3596320_GM18502_barcodes.tsv.gz", col_names = "barcode") %>%
    mutate(barcode = str_remove(barcode, "-1"))
  qc <- read.delim("../data/GSE126321/GSM3596320_GM18502_cellQC.tsv.gz", sep = "\t") %>%
    rownames_to_column("barcode")  %>%
    as_tibble()
  col_df <- left_join(barcodes, qc, by = "barcode") %>%
    as.data.frame()
  gene_info <- AnnotationHub::AnnotationHub()[["AH53537"]] %>%
    as.data.frame() %>%
    group_by(gene_id) %>%
    summarize(chromosome = dplyr::first(seqnames),
              gene_name = dplyr::first(gene_name),
              gene_biotype = dplyr::first(gene_biotype))
  row_df <- genes %>%
    left_join((gene_info), by = "gene_id") %>%
    transmute(gene_id, gene_name = gene_name.x, chromosome, gene_biotype) %>%
    as.data.frame()
  count_mat <- as(mat, "dgCMatrix")
  rownames(count_mat) <- row_df$gene_id
  colnames(count_mat) <- col_df$barcode
  saveRDS(SummarizedExperiment(S4Vectors::SimpleList(counts = count_mat), colData = col_df, rowData = row_df), "../data/GSE126321.Rds")
}

Technical control experiments

se <- zellkonverter::readH5AD("../data/klein_2015.h5ad")
sf <- colSums2(assay(se))
thres <- quantile(sf, 0.5) * c(1, 1.3)
hist(sf, breaks = 50); abline(v = thres, col = "red", lwd = 2)

se_red <- se[, sf > thres[1] & sf < thres[2]]
sf_red <- sf[sf > thres[1] & sf < thres[2]]
se_red <- se_red[rowSums2(assay(se_red)) > 0]
dim(se_red)

[1] 25213   424

mu <- rowMeans2(assay(se_red))
var <- rowVars(assay(se_red))

tibble(mu, var, ercc_gene = str_starts(rownames(se_red), "ERCC-")) %>%
  mutate(mu_pred = mu, mu_sq_pred = mu^2) %>%
  ggplot(aes(x = mu, y = var)) +
    geom_line(data = poisson_pred, aes(group = factor), color = "lightgray", size = 0.1) +
    geom_line(data = gampoi_pred %>% filter((factor == 100 & var < 1.5e3) | (factor != 100 & var < 4e3)), 
              aes(group = factor), color = "#DEB554", size = 0.8) +
    geom_line(data = poisson_pred %>% filter(factor == 1 & mu < 4e3), aes(group = factor), color = "#C981DE", size = 1.2) +
    coord_fixed(expand = FALSE, clip = "off") +
    ggrastr::geom_point_rast(aes(color = ercc_gene), size = 0.1, show.legend = FALSE) +
    annotate(shadowtext:::GeomShadowText, x = 5e3, y = 5e3, label = expression(Var==mu),
             hjust = 0, vjust = 0.5, angle = 45, size = 4, color = "black", bg.colour = "white") +
    annotate(shadowtext:::GeomShadowText, x = sqrt(4e3), y = 4e3, label = expression(Var==mu+1*mu^2),
             hjust = 0, vjust = 0.4, angle = atan(2) / pi * 180, size = 4, color = "black", bg.colour = "white") +
    annotate(shadowtext:::GeomShadowText, x = sqrt(5e3 / 0.01), y = 5e3, label = expression(Var==mu+0.01*mu^2),
             hjust = 0, vjust = 0.4, angle = atan(2) / pi * 180, size = 4, color = "black", bg.colour = "white") +
    annotate(shadowtext:::GeomShadowText, x = sqrt(2e3 / 100), y = 2e3, label = expression(Var==mu+100*mu^2),
             hjust = 0, vjust = 0.4, angle = atan(2) / pi * 180, size = 4, color = "black", bg.colour = "white") +
    annotation_logticks(scaled = TRUE, outside = FALSE, size = 0.2,
                        short = unit(0.05, "cm"), mid = unit(0.1, "cm"), long = unit(0.15, "cm")) +
    scale_x_log10(breaks = c(0.001, 0.1, 10, 1000, 1e5), limits = c(1e-3, 1e5),
                  labels = c(expression(10^-3), expression(10^-1), 
                             expression(10),  expression(10^3), expression(10^5)),
                  name = expression(Mean~(mu))) +
    scale_y_log10(breaks = c(0.001, 0.1, 10, 1000, 1e5), limits = c(1e-3, 1e6),
                  labels = c(expression(10^-3), expression(10^-1), 
                             expression(10),  expression(10^3), expression(10^5)),
                  name = expression(Variance)) +
    scale_color_manual(values = c("TRUE" = "#6c58d1", "FALSE" = "black")) +
    ggtitle("Klein 2015", subtitle = "<span style = 'color: #6c58d1'>ERCC spike-ins</span>") +
    theme(plot.subtitle = ggtext::element_markdown())

Warning: Removed 9803 row(s) containing missing values (geom_path).
Warning: Removed 300 row(s) containing missing values (geom_path).
Warning: Removed 100 row(s) containing missing values (geom_path).
Warning in is.na(x) :
  is.na() applied to non-(list or vector) of type 'expression'
Warning in is.na(x) :
  is.na() applied to non-(list or vector) of type 'expression'
Warning in is.na(x) :
  is.na() applied to non-(list or vector) of type 'expression'
Warning in is.na(x) :
  is.na() applied to non-(list or vector) of type 'expression'

tech_p1 <- last_plot()

se <- zellkonverter::readH5AD("../data/svensson_2017_1.h5ad")
sf <- colSums2(assay(se))
thres <- quantile(sf, 0.5) * c(1, 1.3)
hist(sf, breaks = 50); abline(v = thres, col = "red", lwd = 2)

se_red <- se[, sf > thres[1] & sf < thres[2]]
sf_red <- sf[sf > thres[1] & sf < thres[2]]
se_red <- se_red[rowSums2(assay(se_red)) > 0]
dim(se_red)

[1] 17906   894

mu <- rowMeans2(assay(se_red))
var <- rowVars(assay(se_red))

tibble(mu, var, ercc_gene = str_starts(rownames(se_red), "ERCC-")) %>%
  mutate(mu_pred = mu, mu_sq_pred = mu^2) %>%
  ggplot(aes(x = mu, y = var)) +
    geom_line(data = poisson_pred, aes(group = factor), color = "lightgray", size = 0.1) +
    # geom_line(data = gampoi_pred %>% filter(var < 4e3), aes(group = factor), color = "#DEB554", size = 0.8) +
    # geom_line(data = poisson_pred %>% filter(factor == 1 & mu < 4e3), aes(group = factor), color = "#C981DE", size = 1.2) +
    geom_line(data = gampoi_pred, aes(group = factor), color = "#DEB554", size = 0.8) +
    geom_line(data = poisson_pred %>% filter(factor == 1), aes(group = factor), color = "#C981DE", size = 1.2) +
    coord_fixed(expand = FALSE, clip = "off") +
    ggrastr::geom_point_rast(aes(color = ercc_gene), size = 0.1, show.legend = FALSE) +
    annotation_logticks(scaled = TRUE, outside = FALSE, size = 0.2,
                        short = unit(0.05, "cm"), mid = unit(0.1, "cm"), long = unit(0.15, "cm")) +
    scale_x_log10(breaks = c(0.001, 0.1, 10, 1000, 1e5), limits = c(1e-3, 1e5),
                  labels = c(expression(10^-3), expression(10^-1), 
                             expression(10),  expression(10^3), expression(10^5)),
                  name = expression(Mean~(mu))) +
    scale_y_log10(breaks = c(0.001, 0.1, 10, 1000, 1e5), limits = c(1e-3, 1e6),
                  labels = c(expression(10^-3), expression(10^-1), 
                             expression(10),  expression(10^3), expression(10^5)),
                  name = expression(Variance)) +
    scale_color_manual(values = c("TRUE" = "#6c58d1", "FALSE" = "black")) +
    ggtitle("Svensson 2017 (1)", subtitle = "") +
    theme(plot.subtitle = ggtext::element_markdown())

Warning: Removed 9803 row(s) containing missing values (geom_path).
Warning: Removed 1203 row(s) containing missing values (geom_path).
Warning: Removed 200 row(s) containing missing values (geom_path).

tech_p2 <- last_plot()

se <- zellkonverter::readH5AD("../data/svensson_2017_2.h5ad")
sf <- colSums2(assay(se))
thres <- quantile(sf, 0.5) * c(1, 1.3)
hist(sf, breaks = 50); abline(v = thres, col = "red", lwd = 2)

se_red <- se[, sf > thres[1] & sf < thres[2]]
sf_red <- sf[sf > thres[1] & sf < thres[2]]
se_red <- se_red[rowSums2(assay(se_red)) > 0]
dim(se_red)

[1] 18812   803

mu <- rowMeans2(assay(se_red))
var <- rowVars(assay(se_red))

tibble(mu, var, ercc_gene = str_starts(rownames(se_red), "ERCC-")) %>%
  mutate(mu_pred = mu, mu_sq_pred = mu^2) %>%
  ggplot(aes(x = mu, y = var)) +
    geom_line(data = poisson_pred, aes(group = factor), color = "lightgray", size = 0.1) +
    # geom_line(data = gampoi_pred %>% filter(var < 4e3), aes(group = factor), color = "#DEB554", size = 0.8) +
    # geom_line(data = poisson_pred %>% filter(factor == 1 & mu < 4e3), aes(group = factor), color = "#C981DE", size = 1.2) +
    geom_line(data = gampoi_pred, aes(group = factor), color = "#DEB554", size = 0.8) +
    geom_line(data = poisson_pred %>% filter(factor == 1), aes(group = factor), color = "#C981DE", size = 1.2) +
    coord_fixed(expand = FALSE, clip = "off") +
    ggrastr::geom_point_rast(aes(color = ercc_gene), size = 0.1, show.legend = FALSE) +
    annotation_logticks(scaled = TRUE, outside = FALSE, size = 0.2,
                        short = unit(0.05, "cm"), mid = unit(0.1, "cm"), long = unit(0.15, "cm")) +
    scale_x_log10(breaks = c(0.001, 0.1, 10, 1000, 1e5), limits = c(1e-3, 1e5),
                  labels = c(expression(10^-3), expression(10^-1), 
                             expression(10),  expression(10^3), expression(10^5)),
                  name = expression(Mean~(mu))) +
    scale_y_log10(breaks = c(0.001, 0.1, 10, 1000, 1e5), limits = c(1e-3, 1e6),
                  labels = c(expression(10^-3), expression(10^-1), 
                             expression(10),  expression(10^3), expression(10^5)),
                  name = expression(Variance)) +
    scale_color_manual(values = c("TRUE" = "#6c58d1", "FALSE" = "black")) +
    ggtitle("Svensson 2017 (2)", subtitle = "") +
    theme(plot.subtitle = ggtext::element_markdown())

Warning: Removed 9803 row(s) containing missing values (geom_path).
Warning: Removed 1203 row(s) containing missing values (geom_path).
Warning: Removed 200 row(s) containing missing values (geom_path).

tech_p3 <- last_plot()

Biological control

se <- zellkonverter::readH5AD("../data/nih3t3.h5ad")
sf <- colSums2(assay(se))
thres <- quantile(sf, 0.5) * c(1, 1.3)
hist(sf, breaks = 50); abline(v = thres, col = "red", lwd = 2)

se_red <- se[, sf > thres[1] & sf < thres[2]]
sf_red <- sf[sf > thres[1] & sf < thres[2]]
se_red <- se_red[rowSums2(assay(se_red)) > 0]
dim(se_red)

[1] 19406   788

mu <- rowMeans2(assay(se_red))
var <- rowVars(assay(se_red))

tibble(mu, var, cell_cycle_gene = rownames(se_red) %in% paste0("mm10_", mouse_cell_cycle_genes)) %>%
  mutate(mu_pred = mu, mu_sq_pred = mu^2) %>%
  ggplot(aes(x = mu, y = var)) +
    geom_line(data = poisson_pred, aes(group = factor), color = "lightgray", size = 0.1) +
    geom_line(data = gampoi_pred %>% filter((factor == 100 & var < 1.5e3) | (factor != 100 & var < 4e3)), 
              aes(group = factor), color = "#DEB554", size = 0.8) +
    geom_line(data = poisson_pred %>% filter(factor == 1 & mu < 4e3), aes(group = factor), color = "#C981DE", size = 1.2) +
    # geom_line(data = gampoi_pred, aes(group = factor), color = "#DEB554", size = 0.8) +
    # geom_line(data = poisson_pred %>% filter(factor == 1), aes(group = factor), color = "#C981DE", size = 1.2) +
    coord_fixed(expand = FALSE, clip = "off") +
    # geom_point(size = 0.3, alpha = 0.3) +
    ggrastr::geom_point_rast(aes(color = cell_cycle_gene), size = 0.1, show.legend = FALSE) +
    annotate(shadowtext:::GeomShadowText, x = 5e3, y = 5e3, label = expression(Var==mu),
             hjust = 0, vjust = 0.5, angle = 45, size = 4, color = "black", bg.colour = "white") +
    annotate(shadowtext:::GeomShadowText, x = sqrt(4e3), y = 4e3, label = expression(Var==mu+1*mu^2),
             hjust = 0, vjust = 0.4, angle = atan(2) / pi * 180, size = 4, color = "black", bg.colour = "white") +
    annotate(shadowtext:::GeomShadowText, x = sqrt(5e3 / 0.01), y = 5e3, label = expression(Var==mu+0.01*mu^2),
             hjust = 0, vjust = 0.4, angle = atan(2) / pi * 180, size = 4, color = "black", bg.colour = "white") +
    annotate(shadowtext:::GeomShadowText, x = sqrt(2e3 / 100), y = 2e3, label = expression(Var==mu+100*mu^2),
             hjust = 0, vjust = 0.4, angle = atan(2) / pi * 180, size = 4, color = "black", bg.colour = "white") +
    annotation_logticks(scaled = TRUE, outside = FALSE, size = 0.2,
                        short = unit(0.05, "cm"), mid = unit(0.1, "cm"), long = unit(0.15, "cm")) +
    scale_x_log10(breaks = c(0.001, 0.1, 10, 1000, 1e5), limits = c(1e-3, 1e5),
                  labels = c(expression(10^-3), expression(10^-1), 
                             expression(10),  expression(10^3), expression(10^5)),
                  name = expression(Mean~(mu))) +
    scale_y_log10(breaks = c(0.001, 0.1, 10, 1000, 1e5), limits = c(1e-3, 1e6),
                  labels = c(expression(10^-3), expression(10^-1), 
                             expression(10),  expression(10^3), expression(10^5)),
                  name = expression(Variance)) +
    scale_color_manual(values = c("TRUE" = "#87172f", "FALSE" = "black")) +
    ggtitle("NIH/3T3 Cells", subtitle = "<span style = 'color: #87172f'>Cell cycle marker</span> genes.") +
    theme(plot.subtitle = ggtext::element_markdown())

Warning: Removed 9803 row(s) containing missing values (geom_path).
Warning: Removed 300 row(s) containing missing values (geom_path).
Warning: Removed 100 row(s) containing missing values (geom_path).
Warning in is.na(x) :
  is.na() applied to non-(list or vector) of type 'expression'
Warning in is.na(x) :
  is.na() applied to non-(list or vector) of type 'expression'
Warning in is.na(x) :
  is.na() applied to non-(list or vector) of type 'expression'
Warning in is.na(x) :
  is.na() applied to non-(list or vector) of type 'expression'

p1 <- last_plot()

se <- zellkonverter::readH5AD("../data/hek293t.h5ad")

sf <- colSums2(assay(se))
thres <- quantile(sf, 0.5) * c(1, 1.3)
hist(sf, breaks = 50); abline(v = thres, col = "red", lwd = 2)

se_red <- se[, sf > thres[1] & sf < thres[2]]
sf_red <- sf[sf > thres[1] & sf < thres[2]]
se_red <- se_red[rowSums2(assay(se_red)) > 0]
dim(se_red)

[1] 22804   565

mu <- rowMeans2(assay(se_red))
var <- rowVars(assay(se_red))

tibble(mu, var, cell_cycle_gene = rownames(se_red) %in% paste0("hg19_", human_cell_cycle_genes)) %>%
  mutate(mu_pred = mu, mu_sq_pred = mu^2) %>%
  ggplot(aes(x = mu, y = var)) +
    geom_line(data = poisson_pred, aes(group = factor), color = "lightgray", size = 0.1) +
    # geom_line(data = gampoi_pred %>% filter(var < 4e3), aes(group = factor), color = "#DEB554", size = 0.8) +
    # geom_line(data = poisson_pred %>% filter(factor == 1 & mu < 4e3), aes(group = factor), color = "#C981DE", size = 1.2) +
    geom_line(data = gampoi_pred, aes(group = factor), color = "#DEB554", size = 0.8) +
    geom_line(data = poisson_pred %>% filter(factor == 1), aes(group = factor), color = "#C981DE", size = 1.2) +
    coord_fixed(expand = FALSE, clip = "off") +
    ggrastr::geom_point_rast(aes(color = cell_cycle_gene), size = 0.1, show.legend = FALSE) +
    annotation_logticks(scaled = TRUE, outside = FALSE, size = 0.2,
                        short = unit(0.05, "cm"), mid = unit(0.1, "cm"), long = unit(0.15, "cm")) +
    scale_x_log10(breaks = c(0.001, 0.1, 10, 1000, 1e5), limits = c(1e-3, 1e5),
                  labels = c(expression(10^-3), expression(10^-1), 
                             expression(10),  expression(10^3), expression(10^5)),
                  name = expression(Mean~(mu))) +
    scale_y_log10(breaks = c(0.001, 0.1, 10, 1000, 1e5), limits = c(1e-3, 1e6),
                  labels = c(expression(10^-3), expression(10^-1), 
                             expression(10),  expression(10^3), expression(10^5)),
                  name = expression(Variance)) +
    scale_color_manual(values = c("TRUE" = "#87172f", "FALSE" = "black")) +
    ggtitle("HEK 293T Cells", subtitle = "") +
    theme(plot.subtitle = ggtext::element_markdown())

Warning: Removed 9803 row(s) containing missing values (geom_path).
Warning: Removed 1203 row(s) containing missing values (geom_path).
Warning: Removed 200 row(s) containing missing values (geom_path).

p2 <- last_plot()

se_full <- readRDS("../data/nci-h1975.Rds")
table(se_full$cell_line)

 H1975  H2228 HCC827 
   313    315    274 

se <- se_full[, se_full$cell_line == "H1975"]

sf <- colSums2(assay(se))
thres <- quantile(sf, 0.5) * c(1, 1.3)
hist(sf, breaks = 50); abline(v = thres, col = "red", lwd = 2)

se_red <- se[, sf > thres[1] & sf < thres[2]]
sf_red <- sf[sf > thres[1] & sf < thres[2]]
se_red <- se_red[rowSums2(assay(se_red)) > 0]
dim(se_red)

[1] 15932    99

mu <- rowMeans2(assay(se_red))
var <- rowVars(assay(se_red))

tibble(mu, var, cell_cycle_gene = rownames(se_red) %in% human_cell_cycle_genes) %>%
  mutate(mu_pred = mu, mu_sq_pred = mu^2) %>%
  ggplot(aes(x = mu, y = var)) +
    geom_line(data = poisson_pred, aes(group = factor), color = "lightgray", size = 0.1) +
    # geom_line(data = gampoi_pred %>% filter(var < 4e3), aes(group = factor), color = "#DEB554", size = 0.8) +
    # geom_line(data = poisson_pred %>% filter(factor == 1 & mu < 4e3), aes(group = factor), color = "#C981DE", size = 1.2) +
    geom_line(data = gampoi_pred, aes(group = factor), color = "#DEB554", size = 0.8) +
    geom_line(data = poisson_pred %>% filter(factor == 1), aes(group = factor), color = "#C981DE", size = 1.2) +
    coord_fixed(expand = FALSE, clip = "off") +
    ggrastr::geom_point_rast(aes(color = cell_cycle_gene), size = 0.1, show.legend = FALSE) +
    annotation_logticks(scaled = TRUE, outside = FALSE, size = 0.2,
                        short = unit(0.05, "cm"), mid = unit(0.1, "cm"), long = unit(0.15, "cm")) +
    scale_x_log10(breaks = c(0.001, 0.1, 10, 1000, 1e5), limits = c(1e-3, 1e5),
                  labels = c(expression(10^-3), expression(10^-1), 
                             expression(10),  expression(10^3), expression(10^5)),
                  name = expression(Mean~(mu))) +
    scale_y_log10(breaks = c(0.001, 0.1, 10, 1000, 1e5), limits = c(1e-3, 1e6),
                  labels = c(expression(10^-3), expression(10^-1), 
                             expression(10),  expression(10^3), expression(10^5)),
                  name = expression(Variance)) +
    scale_color_manual(values = c("TRUE" = "#87172f", "FALSE" = "black")) +
    ggtitle("NCI-H1975 Cells", subtitle = "") +
    theme(plot.subtitle = ggtext::element_markdown())

Warning: Removed 9803 row(s) containing missing values (geom_path).
Warning: Removed 1203 row(s) containing missing values (geom_path).
Warning: Removed 200 row(s) containing missing values (geom_path).

p3 <- last_plot()

se <- readRDS("../data/GSE126321.Rds")
# table(se$cellPhase)
# se <- se[, se$cellPhase == "S"]

sf <- colSums2(assay(se))
thres <- quantile(sf, 0.5) * c(1, 1.3)
hist(sf, breaks = 50); abline(v = thres, col = "red", lwd = 2)

se_red <- se[, sf > thres[1] & sf < thres[2]]
sf_red <- sf[sf > thres[1] & sf < thres[2]]
se_red <- se_red[rowSums2(assay(se_red)) > 0]
dim(se_red)

[1] 17909  1242

mu <- rowMeans2(assay(se_red))
var <- rowVars(assay(se_red))

tibble(mu, var, cell_cycle_gene = rownames(se_red) %in% human_cell_cycle_genes) %>%
  mutate(mu_pred = mu, mu_sq_pred = mu^2) %>%
  ggplot(aes(x = mu, y = var)) +
    geom_line(data = poisson_pred, aes(group = factor), color = "lightgray", size = 0.1) +
    # geom_line(data = gampoi_pred %>% filter(var < 4e3), aes(group = factor), color = "#DEB554", size = 0.8) +
    # geom_line(data = poisson_pred %>% filter(factor == 1 & mu < 4e3), aes(group = factor), color = "#C981DE", size = 1.2) +
    geom_line(data = gampoi_pred, aes(group = factor), color = "#DEB554", size = 0.8) +
    geom_line(data = poisson_pred %>% filter(factor == 1), aes(group = factor), color = "#C981DE", size = 1.2) +
    coord_fixed(expand = FALSE, clip = "off") +
    ggrastr::geom_point_rast(aes(color = cell_cycle_gene), size = 0.1, show.legend = FALSE) +
    annotation_logticks(scaled = TRUE, outside = FALSE, size = 0.2,
                        short = unit(0.05, "cm"), mid = unit(0.1, "cm"), long = unit(0.15, "cm")) +
    scale_x_log10(breaks = c(0.001, 0.1, 10, 1000, 1e5), limits = c(1e-3, 1e5),
                  labels = c(expression(10^-3), expression(10^-1), 
                             expression(10),  expression(10^3), expression(10^5)),
                  name = expression(Mean~(mu))) +
    scale_y_log10(breaks = c(0.001, 0.1, 10, 1000, 1e5), limits = c(1e-3, 1e6),
                  labels = c(expression(10^-3), expression(10^-1), 
                             expression(10),  expression(10^3), expression(10^5)),
                  name = expression(Variance)) +
    scale_color_manual(values = c("TRUE" = "#87172f", "FALSE" = "black")) +
    ggtitle("GM18502 Cells", subtitle = "") +
    theme(plot.subtitle = ggtext::element_markdown())

Warning: Removed 9803 row(s) containing missing values (geom_path).
Warning: Removed 1203 row(s) containing missing values (geom_path).
Warning: Removed 200 row(s) containing missing values (geom_path).
Warning: Removed 1812 rows containing missing values (geom_point).

p4 <- last_plot()

p_tech <- cowplot::plot_grid(tech_p1, tech_p2, tech_p3, NULL, nrow = 1, align = "h")

Warning: Removed 9803 row(s) containing missing values (geom_path).
Warning: Removed 300 row(s) containing missing values (geom_path).
Warning: Removed 100 row(s) containing missing values (geom_path).
Warning in is.na(x) :
  is.na() applied to non-(list or vector) of type 'expression'
Warning in is.na(x) :
  is.na() applied to non-(list or vector) of type 'expression'
Warning in is.na(x) :
  is.na() applied to non-(list or vector) of type 'expression'
Warning in is.na(x) :
  is.na() applied to non-(list or vector) of type 'expression'
Warning: Removed 9803 row(s) containing missing values (geom_path).
Warning: Removed 1203 row(s) containing missing values (geom_path).
Warning: Removed 200 row(s) containing missing values (geom_path).
Warning: Removed 9803 row(s) containing missing values (geom_path).
Warning: Removed 1203 row(s) containing missing values (geom_path).
Warning: Removed 200 row(s) containing missing values (geom_path).

p_bio <- cowplot::plot_grid(p1, p2, p3, p4, nrow = 1, align = "h")

Warning: Removed 9803 row(s) containing missing values (geom_path).
Warning: Removed 300 row(s) containing missing values (geom_path).
Warning: Removed 100 row(s) containing missing values (geom_path).
Warning in is.na(x) :
  is.na() applied to non-(list or vector) of type 'expression'
Warning in is.na(x) :
  is.na() applied to non-(list or vector) of type 'expression'
Warning in is.na(x) :
  is.na() applied to non-(list or vector) of type 'expression'
Warning in is.na(x) :
  is.na() applied to non-(list or vector) of type 'expression'
Warning: Removed 9803 row(s) containing missing values (geom_path).
Warning: Removed 1203 row(s) containing missing values (geom_path).
Warning: Removed 200 row(s) containing missing values (geom_path).
Warning: Removed 9803 row(s) containing missing values (geom_path).
Warning: Removed 1203 row(s) containing missing values (geom_path).
Warning: Removed 200 row(s) containing missing values (geom_path).
Warning: Removed 9803 row(s) containing missing values (geom_path).
Warning: Removed 1203 row(s) containing missing values (geom_path).
Warning: Removed 200 row(s) containing missing values (geom_path).
Warning: Removed 1812 rows containing missing values (geom_point).

tech_title <- cowplot::ggdraw() + cowplot::draw_label("(A) Droplets with RNA solution (technical control)", 
                                                      fontface = "bold", x = 0.02, hjust = 0, size = 18)
bio_title <- cowplot::ggdraw() + cowplot::draw_label("(B) Cell line populations (biological control)", 
                                                     fontface = "bold", x = 0.02, hjust = 0, size = 18)

p_res <- cowplot::plot_grid(tech_title, p_tech, bio_title, p_bio, ncol = 1,
                            rel_heights = c(0.2, 1, 0.2, 1))
p_res

cowplot::save_plot("../output/mean_variance_relation_homnogeneous_cells.pdf", p_res, nrow = 2, ncol = 4, base_asp = 0.68, base_height = 5)

Session Info

sessionInfo()

R version 4.1.0 (2021-05-18)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Mojave 10.14.6

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] org.Mm.eg.db_3.13.0         org.Hs.eg.db_3.13.0         AnnotationDbi_1.54.1        SingleCellExperiment_1.14.1 SummarizedExperiment_1.22.0 Biobase_2.52.0             
 [7] GenomicRanges_1.44.0        GenomeInfoDb_1.28.0         IRanges_2.26.0              S4Vectors_0.30.0            BiocGenerics_0.38.0         MatrixGenerics_1.4.0       
[13] matrixStats_0.59.0          forcats_0.5.1               stringr_1.4.0               dplyr_1.0.7                 purrr_0.3.4                 readr_1.4.0                
[19] tidyr_1.1.3                 tibble_3.1.2                ggplot2_3.3.4               tidyverse_1.3.1            

loaded via a namespace (and not attached):
  [1] ggbeeswarm_0.6.0              colorspace_2.0-1              ellipsis_0.3.2                markdown_1.1                  XVector_0.32.0                fs_1.5.0                     
  [7] gridtext_0.1.4                ggtext_0.1.1                  rstudioapi_0.13               farver_2.1.0                  bit64_4.0.5                   interactiveDisplayBase_1.30.0
 [13] fansi_0.5.0                   lubridate_1.7.10              xml2_1.3.2                    sparseMatrixStats_1.4.0       cachem_1.0.5                  knitr_1.33                   
 [19] jsonlite_1.7.2                Cairo_1.5-12.2                broom_0.7.7                   dbplyr_2.1.1                  png_0.1-7                     shiny_1.6.0                  
 [25] BiocManager_1.30.16           compiler_4.1.0                httr_1.4.2                    basilisk_1.4.0                backports_1.2.1               assertthat_0.2.1             
 [31] Matrix_1.3-4                  fastmap_1.1.0                 cli_2.5.0                     later_1.2.0                   htmltools_0.5.1.1             tools_4.1.0                  
 [37] gtable_0.3.0                  glue_1.4.2                    GenomeInfoDbData_1.2.6        rappdirs_0.3.3                Rcpp_1.0.6                    cellranger_1.1.0             
 [43] vctrs_0.3.8                   Biostrings_2.60.1             zellkonverter_1.2.0           xfun_0.24                     rvest_1.0.0                   mime_0.10                    
 [49] lifecycle_1.0.0               AnnotationHub_3.0.0           zlibbioc_1.38.0               scales_1.1.1                  basilisk.utils_1.4.0          ragg_1.1.3                   
 [55] hms_1.1.0                     promises_1.2.0.1              yaml_2.2.1                    curl_4.3.1                    memoise_2.0.0                 reticulate_1.20              
 [61] ggrastr_0.2.3                 stringi_1.6.2                 RSQLite_2.2.7                 BiocVersion_3.13.1            filelock_1.0.2                systemfonts_1.0.2            
 [67] rlang_0.4.11                  pkgconfig_2.0.3               bitops_1.0-7                  lattice_0.20-44               shadowtext_0.0.8              cowplot_1.1.1                
 [73] bit_4.0.4                     tidyselect_1.1.1              magrittr_2.0.1                R6_2.5.0                      generics_0.1.0                DelayedArray_0.18.0          
 [79] DBI_1.1.1                     pillar_1.6.1                  haven_2.4.1                   withr_2.4.2                   KEGGREST_1.32.0               RCurl_1.98-1.3               
 [85] dir.expiry_1.0.0              modelr_0.1.8                  crayon_1.4.1                  utf8_1.2.1                    BiocFileCache_2.0.0           grid_4.1.0                   
 [91] readxl_1.3.1                  blob_1.2.1                    reprex_2.0.0                  digest_0.6.27                 xtable_1.8-4                  httpuv_1.6.1                 
 [97] textshaping_0.3.5             munsell_0.5.0                 beeswarm_0.4.0                vipor_0.4.5